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ABSTRACT 

Near-future cosmological observations targeted at investigations of dark energy pose stringent requirements 
on the accuracy of theoretical predictions for the clustering of matter. Currently, A-body simulations comprise 
the only viable approach to this problem. In this paper we demonstrate that A-body simulations can indeed be 
sufficiently controlled to fulfill these requirements for the needs of ongoing and near-future weak lensing surveys. 
By performing a large suite of cosmological simulation comparison and convergence tests we show that results for 
the nonlinear matter power spectrum can be obtained at 1% accuracy out to k ~ 1 AMpc . The key components 
of these high accuracy simulations are: precise initial conditions, very large simulation volumes, sufficient mass 
resolution, and accurate time stepping. This paper is the first in a series of three, with the final aim to provide a 
high-accuracy prediction scheme for the nonlinear matter power spectrum. 

Subject headings: methods: A-body simulations — cosmology: large-scale structure of universe 



1. INTRODUCTION 

The nature of the dark energy believed to be causing the cur- 
rent accelerated expansion of the Universe is one of the greatest 
puzzles in the physical sciences, with deep implications for our 
understanding of the Universe and fundamental physics. The 
twin aims of better characterizing and further understanding 
the nature of dark energy are widely recognized as key science 
goals for the next decade. Although dark energy remains very 
poorly understood, theory nevertheless plays an essential role 
in furthering this enterprise. 

The phenomenology of cosmological models is theory-driven 
not only in terms of providing explanations for the diverse phe- 
nomena that are observed, as well as promoting alternative ex- 
planations of existing measurements, but also due to the in- 
creasing reliance on theorists to produce sophisticated numer- 
ical models of the Universe which can be used to refine and 
calibrate experimental probes. Without a dedicated effort to de- 
velop the tools and skill-sets necessary for the interpretation of 
the next generation of experiments, we risk being "theory lim- 
ited" in essentially all areas of dark energy studies. 

Forecasts for determination of the dark energy equation of 
state and other cosmological parameters from next-generation 
observations of cosmological structure typically assume cali- 
bration against simulations accurate to the level of 1 % or better. 
This target has rarely been met for simulations of complex non- 
linear phenomena such as the formation of large-scale structure 
in the Universe. However it is precisely these probes, which 
provide information on both the geometry of space-time and 
the growth of large-scale structure, which will be key to unrav- 
eling the mystery of dark energy. 

For upcoming measurements to be exploited to the full, the- 
ory must reach not only the levels of accuracy justified by the 
measurements but also cover a sufficiently wide range of cos- 
mologies. The problem breaks down to two questions: (i) What 
is a reasonable coverage of cosmological parameters, given the 
expected set of observations? (ii) What is the required accuracy 



for theoretical predictions - over this range of parameters - for 
the given set of observations? We fully expect that the answers 
to both (i) and (ii) will evolve, requiring more accurate mod- 
eling of a smaller range of models, so we are most interested 
here in the near-term needs. Associated with the first problem 
is the fact that, given the impossibility of running complex sim- 
ulations over the many thousands of cosmologies necessary for 
grid based or Markov Chain Monte Carlo (MCMC) estimation 
of cosmological parameters, one must develop efficient interpo- 
lation methods for theoretical predictions. These methods must 
of course also satisfy the accuracy requirements of question (ii). 

The control of errors in the underlying t heory for the CMB 
is adequate to an alyze results from Planck (Seljak et al. 2003; 
IWong et alj2008b . This is, however, not the case for predictions 
of gravitational clustering in the nonlinear regime, as is required 
for cluster counts, redshift space distortions, baryon acoustic 
oscillations (BAO) and weak lensing (WL) observations. In 
the case of BAO, the galaxy power spectrum in the quasi-linear 
regime should be known to sub-percent accuracy, and for WL 
the same is true for the mass power spectrum to significantly 
smaller scales. Perturbation theory has errors on the mass 
power spectrum currently estimated to be at the percent level in 
the weakly nonlinea r regime (see, e.g. Jjeong & Komatsul2 006; 
Carlson et a ii 120091 for recent treatments or Bernardeau et al. 



2002 for an earlier review). To reduce these errors, test the 



approximations, and model galaxy bias, numerical simulations 
are unavoidable. Theoretical templates, in terms of current 
power spectrum fits based on simulations (with errors at the 
5% level), are already a limiting factor for WL observations at 
wavenumbers k ~ 1/iMpc -1 . iHuterer & Takadal d2005l) show 
that in order to avoid errors from imprecise theoretical tem- 
plates mimicking the effect of cosmological parameter varia- 
tions, the power spectrum has to be calibrated at about 0.5-1% 
for 0.1 /zMpc -1 <k< 10/iMpc -1 . The scale most sensitive for 
WL measurements is around k ~ 1 /iMpc" 1 and z ~ 0.5 and the 
power spectrum therefore needs to be calibrated the most ac- 
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curately at that point (see, e.g, iHuterer & Takadal (|2005 ), Fig- 
ure 1). In a very recent paper. iHilbert et al] (12008) re-emphasize 
the need for very accurate predictions for the theoretical power 
spectr um, pointing out that curre ntly used fitting functions such 
as the | Peacock & Doddsl (fl996) formula or the fit derived by 
ISmith et all d2003l) underestimate the cosmic shear-power spec- 
tra by > 30% for k > 10 h Mpc -1 . 

In order to extract precise cosmological information from 
WL measurements, additional physics beyond the gravitational 
contribution must be taken into account. At length scales 
smaller than k ~ 1/tMpc , baryo nic effects are expected to 
enter at a competing level ([White! £0041 IZhan & Kn ox 2004; 
iJing etal.ll2006l iRudd et al .1120081) and will have to be treated 
separately, either directly via hydrodynamic simulations, or 
as is more likely, by a combination of simulations and self- 
calibration techniques. In any case, gravitational A^-body sim- 
ulations must remain the bedrock on which all of these tech- 
niques are based. 

Given the success of the CDM paradigm in explaining cur- 
rent observational data we shall consider cosmologies within 
that framework. All of our models will assume a spatially flat 
Universe with purely adiabatic fluctuations and a power-law 
power spectrum. Since it is unlikely that near-term observa- 
tions can place meaningful constraints on the temporal varia- 
tion of the equation of state of the dark energy, we will restrict 
attention to cosmologies with a constant equation of state pa- 
rameter w = -p/p (where p is the pressure and p the density 
of the dark energy with w = -1 in a ACDM cosmology). Since 
ACDM is a good fit to the data, the accuracy of simulations can 
be established primarily around this point. 

Taking all these considerations into account, the purpose of 
this paper is to establish that gravitational A^-body simulations 
can produce P(k) results accurate to 1 % out to k ~ 1 h Mpc -1 
between z = — 1 for dark energy models with a constant w. In 
this regime, additional physics is controllable at the required 
level of accuracy. The target regime covers the most important 
range for current and near-future WL surveys. Showing that the 
required accuracy can be obtained from A^-body simulations is 
only the first step in setting up a power spectrum determination 
scheme useful for weak lensing surveys. In order to analyze 
observational data and infer cosmological parameters, precise 
predictions for the power spectrum over a large range of cos- 
mologies are required. This paper - establishing that achieving 
the base accuracy is possible - is the first in a series of three 
communications. In the second, we will demonstrate that a rel- 
atively small number of numerically obtained power spectra are 
sufficient to derive an accurate prediction scheme - or emula- 
tor - for the power spectrum covering the full range of desired 
cosmologies. The third paper of the series will present results 
from the complete simulation suite, named the "Coyote Uni- 
verse" after the computing cluster on which it has been carried 
out. The third paper will also contain a public release of a pre- 
cision power spectrum emulator. 

In order to establish the accuracy over the required spatial 
dynamic range, as well as over the redshifts probed, a va- 
riety of tests need to be conducted. These include studies 
of the initial conditions, convergence to linear theory at very 
large length scales, the mass resolution requirement, and other 
evolution-specific requirements such as force resolution and 
time-stepping errors. To establish robustness of the final re- 
sults, codes based on different A^-body algorithms should inde- 
pendently converge to the same results (within error bounds). 



While some of these studies have been conducted separately 
and within the confine s of the cosmic code verification project 
dHeitmann et alj 120071) . this is the first time that the more or 
less complete set of possible problems has been investigated in 
realistic simulations. 

We find that it is indeed possible to control the accuracy of 
A^-body simulations at 1% out to k ~ 1 /zMpc" 1 . Even though 
these scales are not very small, the simulation requirements 
are rather demanding. First, the simulation volume needs to 
be large enough to capture the linear regime accurately. Due 
to mode-mode coupling, nonlinear effects influence scales as 
large as 500 /i _I Mpc. Therefore, the simulation volume needs 
to cover at least 1 (/r'Gpc) 3 . Second, with this requirement im- 
posed, the number of particles necessary to avoid errors from 
discreteness effects at the smallest length scales of interest, 
also becomes substantial. As we discuss later, numerical re- 
sults aiming for accuracy at the sub-percent level, can only be 
trust ed at scales below the particle Nyquist wavenumber (see 
also iJovce eTail 120081) . We find that a 1 (/r'Gpc) 3 simula- 
tion volume leads to a minimum particle loading of a billion 
particles. Third, it is important to start the simulation at a 
high enough redshift to allow enough dynamic range (in time) 
for structures to evolve correctly and for the initial perturba- 
tions to be captured accurately by the Zel'dovich approxima- 
tion. Lastly, the force resolution and time stepping has to be ac- 
curate enough to ensure convergence of the simulation results. 

The paper is organized as follows. In Section|2]we use a sim- 
ple example to demonstrate the need for precision predictions 
from theory. Section [3] contains a description of the A^-body 
codes used in this paper and some basic information about the 
simulations. In Section [4] we briefly describe the power spec- 
trum estimator. In Sections [5] and [6] investigations of initial 
conditions and time evolution are reported, demonstrating that 
the simulations can achieve the required accuracy levels. Fi- 
nally, we compare the numerical re sults to the c ommonly used 
semi-analytic HaloFit approach (ISmith et al.l 120031) in Sec- 
tion [7] finding a discrepancy of ~ 5- 10% between the fit and 
the simulations. We provide a summary discussion of our re- 
sults in Section [8] Appendix [A] discusses errors in setting up 
the initial conditions, comparing the Zel'dovich and 2LPT ap- 
proximations. Appendix IB1 provides details of the Richardson 
extrapolation procedure used for some of the convergence tests. 

2. THE PRECISION COSMOLOGY CHALLENGE 

Before discussing how to achieve 1 % accuracy for the non- 
linear power spectrum, we will briefly demonstrate the impor- 
tance of accurately determining the power spectrum. In our 
example, we assume the ability to measure the power spec- 
trum from observations at 1 % accuracy in the quasi-linear and 
nonlinear regimes. On larger scales, accounting for sample 
variance (statistical limitations due to finite volume-sampling) 
leads to an increase in the statistical error, of up to 10%. These 
values are rough estimates, which are sufficient to make our 
point in this simple example. 

We generate two sets of mock measurements: one from a 
power spectrum generated with a halo model-inspired fitting 
formula given by the code HaloFit as implemented in CAMB 
l |http : //camb . info) and another directly from a set of 
high-precision simulations. We then move points off the base 
power spectrum according to a Gaussian distribution with vari- 
ance specified by the error estimates given above. The resulting 
mock data points and the underlying power spectra are shown 
in Figure Q] On a logarithmic scale, the data points and power 
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FIG. 1 . — Upper panel: Synthetic data from a HALOFlT run. Lower panel: 
Synthetic data from a combination of several /V-body runs. In both cases the 
black line shows the underlying power spectrum from which the data was 
drawn and the red points show 34 data points with error bars. At small spatial 
scales, the assumed error is 1%, rising to 10% at large scales due to increased 
sample variance. 



spectra are almost indistinguishable. As we will show later in 
Section |7J the difference between the HaloFit and A^-body 
power spectra is at the 5-10% level: this difference is enough to 
lead to significant biases in parameter estimation. 

We determine the best-fit parameters from the two mock data 
sets using the following parameter priors: 

0.02 < uj b < 0.025, 

0.11 <cj„, <0.15, 

0.85 < n s < 1.05, 

-1.3 < w <-0.7, 

0.7 < a s <0.9, (1) 

where ui/, = VLi,h 2 and u> m = Q m h 2 - We do not treat h as an 
independent variable but determine it via the CMB constraint 
I a = iTdi ss lr s = 302.4 where di ss is the distance to the last scat- 
tering surface and r s is the sound horizon (more details of how 
we construct our model grid will be given in Paper II). 

The parameter estimation analysis then proceeds via a com- 
bination of model interpolation and Markov Chain Monte Carlo 
(MCMC) as implement ed in our recentl y introduced cosmic 
calibration framework dHabibetalJ 120071) . We use HaloFit 
to generate the nonlinear power spectra for the MCMC analy- 
sis. That is, we analyze a HaloFit synthetic data set and one 
generated from numerical simulations against a set of model 
predictions from HaloFit generated power spectra. The re- 
sults, which are all obtained from data at z = 0, are shown in 
Figure [2] The upper panel shows the results from the anal- 
ysis of the HaloFit synthetic data, where the parameter es- 
timation works extremely well, being essentially a consistency 
check for the statistical framework. The result also points to the 
constraining power of matter power spectrum data. The lower 
panel in Figure [2] shows the corresponding result for the syn- 
thetic data generated directly from the simulations. In this case, 
the ~ 5% errors in the HaloFit model predictions are clearly 
seen to be problematic: most of the parameters are significantly 
off, u> m and w being mis-estimated by ~ 20%. 

To illuminate this result further, Figure [3] shows how the dif- 
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FIG. 2. — Posterior distributions for the five parameters under consideration. 
Upper panel: Results for the analysis of the HALOFIT synthetic data set an- 
alyzed with a set of HaloFit power spectra. The red dots indicate the true 
values. As is to be expected, the constraints on the parameters are very good. 
Lower panel: Results for the HALOFlT-based analysis of the W-body synthetic 
data set. Note that the constraints for u) m and w are now incorrect at ~ 20%. 
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FIG. 3. — Influence of the different cosmological parameters on the matter 
power spectrum. The y-axis shows the deviation of the logarithm of the power 
spectrum from its value when each parameter is set at the midpoint of its prior 
range, as specified in Eqn. (T}. In each of the five sub-figures, only one pa- 
rameter is varied between the maximum and the minimum of the prior range, 
while the other four parameters are kept at the mid-point value. The light to 
dark lines correspond to the smallest parameter setting to the largest. Note that 
the Hubble parameter h is not a free parameter in our study but determined 
for each cosmology from CMB constraints. Consequently, k is measured in 
Mpc~' in this plot. The variations of h influence the sensitivity plot for the 
dark energy equation of state, w, most strongly. 



ferent cosmological parameters affect the matter power spec- 
trum. This analysis is based on HaloFit power spectra. (For 
an earlier investigation on the sensitivity of the power spec- 
trum with respect to th e dark energy equation of state, see, e.g., 
iMcDonald et al.ll2006h . In each of the five sub-panels, only one 
parameter is varied while the others are held fixed at the mid- 
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point of their prior range, given in Eqs. ([l} (the value for the 
Hubble parameter h is different for each cosmology and set by 
CMB constraints). The y-axis shows the deviation of the natu- 
ral logarithm from the nominal spectrum. As we show in Sec- 
tion [7] HaloFit underpredicts the nonlinear power spectrum 
even on quasi-linear scales. This can lead to a suppression in 
the prediction for the spectral index n s , inducing in turn an over- 
prediction of u> m , since both parameters have similar influence 
on the power spectrum on large scales. Figure|3]also shows that 
the change due to a% and w is similar at similar scales, in rough 
agreement with the overestimation of <rg and the underestima- 
tion of w. Of course the precise interactions of the parameters 
are more complicated and an incorrect posterior distribution of 
one of the parameters will shift all of the others. 

The example used here is certainly too simplified, relying 
only on large scale structure "observations" and making no at- 
tempt to take into account covariance, degeneracies, other ob- 
servations, etc. For example, including a second observational 
probe such as the cosmic microwave background would pro- 
vide a tighter constraint on er^, reducing the 20% shift in w. 
Nevertheless, the example clearly illustrates the general point 
that to perform an unbiased data analysis the theory underlying 
the analysis framework must match or preferably exceed the 
accuracy of the data. 

3. /V-BODY CODES AND SIMULATIONS 

The numerical computations carried out and analyzed in this 
paper are A^-body simulations that model structure formation 
in an expanding universe assuming that gravity dominates all 
other forces. The phase space density field is sampled by 
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FIG. 4. — Upper panel: Comparison of dimensi onless power spectra f rom 
a handful of different codes, taken from the data of Heitmann et al. 1 200/3) for 
the "LCDMb" box: a ACDM model with Cl m = 0.314, h = 0.71, n s = 0.99, and 
L box = 256 /T 1 Mpc. The force resolution of the two PM codes, MC 2 and PMM, 
and the peak resolution of the FLASH code are roughly a factor of ten lower 
than for the other codes (the different force kernels make a precise comparison 
of the force resolution difficult). The dotted lines show the 1% limit. The high 
force-resolution codes agree to 0(1%) up to k ~ l/iMpc~' despite different 
choices for the force softening and other numerical parameters. Lower panel: 
Comparison of GADGET-2 and ART for a simulation with 1024 3 particles 
and L|j OX = 1 /i~'Gpc. The cosmology is very similar as for our major runs, the 
main difference being the starting redshift which is Zm = 65.66. The agreement 
of the two codes is better than 1% over all scales. 



finite-mass particles and these particles are evolved using self- 
consistent force evaluations. Although the effects of baryons 
and neutrinos are taken into account while setting up initial con- 
ditions, only their gravitational contribution to the ensuing non- 
linear dynamics of structure formation is kept (along with that 
of the dark matter). Gas dynamics, feedback effects, etc. are all 
neglected. At sufficiently small scales this neglect is clearly not 
justified, but at the 1% level and for wavenumbers smaller than 
k ~ 1 hMpc~ l this assumption is expected to hold. 

In order to solve the A^-body problem, we employ two com- 
monly used algorithms, the particle-mesh (PM) approach and 
the tree-PM approach. The A^-body methods model many-body 
evolution problems by solving the equations of motion of a 
set of tracer particles which represent a sampling of the sys- 
tem phase space distribution. In PM codes, a computational 
grid is used to increase the efficiency of the self-consistent 
inter-particle force calculation. In the codes used in this pa- 
per, the Vlasov-Poisson system of equations for an expand- 
ing universe is solved using Cloud-in-Cell (CIC) mass depo- 
sition and interpolation with second order (global) symplec- 
tic time-stepping and a Fast Fourier Transform (FFT)-based 
Poisson solver. The advantage of the PM method is good er- 
ror control and speed, the major disadvantage is the restriction 
on force resolution imposed by the biggest FFT that can be 
performed (typical current limits being 2048 3 grids or 4096 3 
grids). Two independently written PM codes were checked 
against each othe r in the low k regime, o ne being the PM code 
MC 2 described in Heit mann et al. (12005b . with excellent agree- 
ment being achieved. In a ddition, the publicly available code 
GADGET-2 (ISpringell2005l) was slightly modified to run in pure 
PM mode. The agreement between these codes was excellent. 

Tree-PM is a hybrid algorithm that combines a long-range 
force computation using a grid-based technique, with shorter- 
range force computation handled by a tree algorithm. The tree 
algorithm is based on the idea that the gravitational potential of 
a far-away group of particles is accurately given by a low-order 
multipole expansion. Particles are first arranged in a hierarchi- 
cal system of groups in a tree structure. Computing the poten- 
tial at a point turns into a descent through the tree. For most 
of our high-resolution runs we use the tree-PM code GADGET- 
2, for some of the tests and c omparison w e also use the code 
TreePM which is described in Whiji (|2002|) . 

Several di fferent A/-body codes have b een compared in pre- 
vious work (iHeitmann et alj|2005l 120071) . including PM, tree- 
PM, adaptive-mesh-refinement, pure tree, and particle-particle 
PM codes. The results of these code verification tests are con- 
sistent with the idea that 1% error control is possible up to 
k ~ 1 /iMpc" 1 (at z = 0), as shown in Figure|4] The upper panel 
in the figure shows a compar ison of the power spe ctra from 
a subset of the codes used in IHeitmann et alj (2007) with re- 
spect to a GADGET-2 run. The simulations are performed with 
256 3 particles in a 256 /r'Mpc box. We find agreement at the 
one-percent level between the high resolution codes despite the 
use of different choices for the force softening and other nu- 
merical parameters. In a separate test, we comp ared GADGET-2 
with the Adaptive Refinement Tr ee (ART) code ( Kra vtsov et alj 
ll997tlGotflober & K lypin 2008). The simulation encompassed 
a volume of (1 /r'Gpc) 3 and 1024 3 particles. The agreement 
between the two codes was again better than one percent be- 
tween z = and z = 1 and out to k ~ 1 /zMpc" 1 . The result for 
Z = is shown in the lower panel of Figure [4] The excellent 
and robust - w.r.t. numerical parameter choices - agreement 
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between different codes provides confidence that it is possible 
to predict the matter power spectrum at the desired accuracy. 

We use a combination of PM and tree-PM runs for this pa- 
per, and in the follow-up work, to create an accurate prediction 
for the matter power spectrum. At quasilinear spatial scales - 
large, yet not fully described by linear theory (k ~ 0.1 /zMpc" 1 ) 
- lower resolution PM simulations are adequate. Furthermore, 
to reduce the variance due to finite volume-sampling - a prob- 
lem at low values of k - simulations should be run with many 
realizations of the same cosmology. We fulfill this requirement 
by running a large number of PM simulations with either 512 3 
or 1024 3 particles. In order to resolve the high-£ part of the 
power spectrum, we use the GADGET-2 code. 

The codes are run with different settings as explicitly dis- 
cussed in the tests mentioned below. In the case of the GADGET- 
2 runs, we use a PM grid twice as large, in each dimension, 
as the number of particles, and a (Gaussian) smoothing of 1 .5 
grid cells. The force matching is set to 6 times the smoothing 
scale, the tree opening criterion being set to 0.5%. The soften- 
ing length is set to 50 kpc. For more general details on t he cod e 
settings in GADGET-2 and the code itself, see Springel (120051) . 

The pure PM simulations have twice as many mesh points in 
each dimension as there are particles. The integration variables 
are the position and conjugate momentum, with time-stepping 
being in constant steps of A In a = 0.02. The forces are ob- 
tained using 4 th order differencing from a potential field com- 
puted using Fourier transforms. The input density field is ob- 
tained frojnJhe_r2articJe_d2Stribution using CIC charge assign- 
ment dHockne v & East wood! 1 19891) and the potential is com- 
puted using a 1 /k 2 kernel. 

If not stated otherwise, our fiducial ACDM model has the 
following cosmological parameters: f2 m = 0.25 for the total 
matter content, a cosmological constant contribution specified 
by tt\ = 0.75, baryon density as set by u>i, = ttbh 2 = 0.024, a 
dimensionless Hubble constant of h = 0.72, the normalization 
specified by a% = 0.8, and a fixed spectral index, n s = 0.97. 
These parameters are in accord with the latest WMAP re- 
sults (Dunkle v et al.l 12008b . The model is run with box size 
of (936 /r'Mpc) 3 and with 1024 3 particles. For some of the 
tests we use a downscaled version of this simulation but keep 
the interparticle spacing approximately the same (1 /r'Mpc). 

4. POWER SPECTRUM ESTIMATION 

The key statistical observable in this paper is the density fluc- 
tuation power spectrum P{k), the Fourier transform of the two- 
point density correlation function. In dimensionless form, the 
power spectrum may be written as 



A 2 (k) : 



k 3 P(k) 



(2) 



which is the contribution to the variance of the density pertur- 
bations per Ink. 

Because TY-body simulations use particles, one does not di- 
rectly compute P(k) or equivalently, A 2 (k). Our procedure is 
to first define a density field on a grid with a fine enough reso- 
lution such that the grid filtering scale is much higher than the 
k scale of interest. This particle deposition step is carried out 
using CIC assignment. The application of a discrete Fourier 
transform (FFT) then yields S(k) from which we can compute 
f(k) = |5(k)| 2 , which in turn can be binned in amplitudes to fi- 
nally obtain P(k). Since the CIC assignment scheme is in effect 
a spatial filter, the smoothing can be compensated by dividing 



P(k) by W 2 (k), where 

w(k)=fJ kxLg 



•2 ( kyLg \ 2 / kyLg 
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(3) 



and L g is the size of the grid cell. Typically the effect of this cor- 
rection is only felt close to the maximum (Nyquist) wavenum- 
ber for the corresponding choice of grid size. One should also 
keep in mind that particle noise and aliasing artifacts can arise 
due to the finite number of particles used in /Y-body simulations 
and due to the finite grid size which is used for the power spec- 
trum estimation. As explained further below, convergence tests 
based on varying the number of sampling particles can help es- 
tablish the smallest length scales at which accurate results can 
be obtained. The particle loading in our simulations is sufficient 
to resolve the power spectrum at the scales of interest, such that 
possible shot noise is at the sub-percent level. 

It is common to make a correction for finite particle number 
by subtracting a Poisson "shot-noise" component from the bin- 
corrected power spectrum: 

3 



Ashot(*)=„ 2 



2vr 2 



N„ 



(4) 



where N p is the cube-root of the number of particles and L is 
the box length. We have not done this in this paper because our 
particle loading is large enough to render it a small correction 
on the scales of interest and it is not clear that this form cap- 
tures the nature of the correction correctly. Note that the initial 
conditions have essentially no shot noise at all, and the evolu- 
tion prior to shell-crossing does not add any. Shot noise thus 
enters through the high-£ sector and filters back to lower A: in a 
complex manner. 

We average P(k) in bins linearly spaced in k of width Ak ~ 
0.001 Mpc -1 , and report this average for each bin containing at 
least one grid point. We assign to each bin the k associated with 
the unweighted average of the k's for each grid point in the bin. 
Note that this procedure introduces a bias in principle, since 
for nonlinear functions (f(x)) ^ f((x)), but our bins are small 
enough to render this bias negligible. 

In a recent paper, IColombi et al.l (120081) suggest a differ- 
ent way to accurately estimate the power spectra from TY-body 
simulations. Their method is based on a Taylor expansion of 
trigonometric functions to replace large FFTs. This way, they 
are able to estimate the power spectrum out to small scales 
with minimal memory overhead, the major obstacle for large 
FFTs. We have checked their method up to fifth order against 
our results from the 2048 3 FFT and found excellent agree- 
ment. Our FFT is clearly large enough to avoid any aliasing 
at ~ 1 /zMpc" 1 . 

5. INITIAL CONDITIONS 

The initial conditions in TY-body codes are often a source of 
systematic error in ways that can sometimes be hard to detect. It 
is, therefore, essential to ensure that the implementation of the 
initial conditions is not a limiting factor in attaining the required 
accuracy of the power spectrum over the redshift range of inter- 
est. An important aspect here is the choice of starting redshift. 
There are two reasons for this: (i) The Lagrangian perturbation 
theory used to generate the initial particle distribution (usually 
the leading order Zel'dovich approximation) is more accurate 
at higher redshifts, and (ii) for a given (nonlinear) k scale of in- 
terest, enough time must have elapsed for the correct nonlinear 
power spectrum to be established at that scale, at the redshift of 
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FIG. 5. — Top two panels: comparison of outputs at z = 10 using different starting redshifts. The particles are colored with respect to their velocities. The 
simulation box is 8 /r'Mpc on a side. The simulation shown in the left panel was started at Zj„ = 250, and the one in the right panel was started at Zm = 50. In the 
simulation that was started z: m = 50, the structures that have formed by z. = 10 are not as concentrated as in simulations with a high-z start, leading to the possible 
lowering of halo masses. The lower panel shows differences along a filament. In this case a line was drawn betwe en each particle pos it ion in the two different dat a 
sets. The longer t he line, the larger the difference due to the two different initial redshifts. For more details see lHaroz et alj (2008); H aroz &Heitmannl (2008): 
ILukic et al.1 f2007l) . 



interest. Note that this has nothing to do w ith the accuracy o f 
the Zel'dovich or some other approximation (ILukic et al.l2 007). 

Due to a combination of the two effects mentioned above, 
delayed starts typically lead to a suppression of structure for- 
mation (including the halo mass function) as shown in Figure|5] 
We now describe our basic methodology for generating initial 
conditions and choosing the starting redshift. 

5.1. Initial Condition Generation 

As is standard, we generate our initial conditions by dis- 
placing particles from a regular Cart esian grid ("quie t start") 
using the Zel'dovich approximation dZerdovichll 19701) . In the 
Zel'dovich approximation, the particle displacement and veloc- 
ity are given by 

x(q) = q-D 1 V^ (I) , (5) 
dn 

y= — = -D l f l HV q ^ l \ (6) 

Here q is the initial (on-grid) position of the particle, x is the 
final position, D\ is the linear growth factor defined below in 
Eqn. ([8]) and (jp^ is the potential field. H is the Hubble con- 
stant, fi is the logarithmic derivative of the growth function 
fi = (dlnDj)/{dlna), and the time-independent potential <fp^ 
obeys the Poisson equation S/ 2 q <fP\q) = S(q). 

A recent suggestion is to determine the initial displacement 
of the particles and their velocities via second order Lagrangian 
perturbation theory (2LPT) instead of using the (leading order) 



Zel'dovich approximation dScoccimarr o 1998L lCrocce et"ail 2006). 
In principle, this could allow a later start of the simulation 
(lower z m ) without losing accuracy in the final result. However, 
it does not address the problem of keeping a sufficient num- 
ber of expansion factors between the initial and final redshifts. 
Additionally, error control of the perturbation theory and its 
convergence properties need to be carefully checked. We have 
therefore decided on a more conservative approach: instead of 
using higher order schemes to generate initial conditions, we 
choose a high enough starting redshift that higher order effects 
are clearly negligible. Since most of the code's runtime is at 
low redshift, the additional overhead for starting the simulation 
early is minimal. Nevertheless, it is informative to compare the 
initial conditions from the 2LPT formalism and the Zel'dovich 
approximation to confirm that the next-order perturbations are 
indeed small. For a detailed comparison of the two approaches, 
see AppendixlAl 

The potential field is generated from a realization of a Gaus- 
sian random density field <5(k) (with random phases). The initial 
power spectrum is 

P(k) = Bk' , T 2 (k), (7) 

where B determines the normalization and T 2 (k) is the matter 
transfer function. We compute T 2 (k) using the numerical code 
CAMB. The results from CAMB were compare d against those 
genera t ed by an independen t co de described in IWhite & Scottl 
(119961) . iHu & Whitd (119971) and iHu et afl (Il998h . The results 
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from this code are known to agree well with CMBfast ( Selia k et ail 
120031) . The final level of agreement was at the ~ 10 3 level for 
the k modes of interest. 

The displacement field is easily generated in Fourier space: 
The Fourier transform of the displacement field is proportional 
to (k/fc 2 )<5(k) in the continuum, and we compute the displace- 
ments using FFTs. The FFT grid is chosen to have twice as 
many points, in each dimension, as there are particles. 

The scale-in dependent linear growth factor, D\(z), satisfies 
dPeacocklll999h 



4 + I^ 

2 p c 



D\ 



3Pm_ £ Pc_ 3 

2 p c 2 p c 



= 0. 

(8) 



Here p c oc H 2 is the critical density, p m is the matter density and 
primes denote differentiation with respect to lna. Our conven- 
tion has D\{z = 0) = 1 and D\(z) oc (1 +z)~ l when p„, ~ p c . This 
procedure neglects the differential evolution of the baryons and 
dark matter, but since we are simulating only collisionless sys- 
tems here this is the most appropriate choice. 

5.2. The Initial Redshift 

The choice of the starting redshift depends on three factors: 
the simulation box size, the particle loading, and the first red- 
shift at which results are desired. The smaller the box and the 
higher the first redshift of interest, the higher the initial red- 
shift must be. It is not easy to provide a universal "recipe" for 
determining the optimal starting redshift. For each simulation 
set-up, convergence tests must be performed for the quantities 
of interest. Nevertheless, there are several guiding principles to 
determine the starting redshift for a given problem. These are: 

• Ensure that any unphysical transients from the initial 
conditions are negligible at the redshift of interest. 

• Ensure a sufficient number of expansion factors to allow 
structures to form correctly at the scales of interest. 

• Ensure that the initial particle move on average is much 
smaller than the initial inter-particle spacing. 

• Ensure that A 2 (k) -^i 1 at the wavenumber of interest. 

A more detailed description - from a mas s funct ion-centric 
point of view - can be found in iLukic et al.l (12007). The aim 
here is to measure the power spectrum from a (936 /r'Mpc) 3 
box between z = 1 and z = at k = 1 /iMpc" 1 at 1 % level accu- 
racy. In order to fulfill the first and second criteria given above, 
we generate the initial conditions at z- m such that D(zm) /D(z = 
1) = 0.01. With Di(z) ~ a(z) = 1/(1 +z) this leads to a starting 
redshift of approximately z; n = 200 and one hundred expansion 
factors between the starting redshift and z = 1 . Note that this 
criterion is completely independent of the box size and particle 
loading, though it is cosmology dependent via the growth rate. 

For the (936/T 1 Mpc) 3 boxes we simulate, this starting red- 
shift leads to rms displacements between 3-5% of the mean 
interparticle spacing, satisfying the condition that the rms dis- 
placement should be much less than the mean interparticle 
spacing. This measurement clearly depends on the box size. 
A smaller box would have led to much bigger displacements 
with respect to the mean inter-particle spacing. At z\ n = 200 the 
dimensionless power at the fundamental mode is O(10~ 8 ) and 
at the Nyquist frequency is 0(1 0~ 4 ) which clearly satisfies the 
last point of the list above. 
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FIG. 6. — Comparison of ratios of the dimensionless power spectra at z = 1 
(upper panel) and z = (lower panel) when evolved using a PM code from 
initial conditions generated using the Zel'dovich approximation at the starting 
redshifts indicated. The rms displacement for the starts is 0.335, 0.168, 0.084, 
and 0.055 times the mean inter-particle spacing (for Zm = 52, 105, 211 and 
317). The dotted lines mark the 1% limit. If the code is started at z\ n — 52, we 
see a suppression of the power spectrum by ~ 3% at z = 1 and ~ 2% at z = 0. 



Finally, in order to verify that the above criteria are sufficient 
to guarantee one percent accuracy between z = 1 and z = we 
carry out a convergence study. As Figure [6] shows, the power 
spectrum is well converged by z = given z m satisfying our cri- 
teria. Our results a re in very g ood agreement with similar tests 
carried out by, e.g. jMal (120071) . We carried out numerous other 
tests with very similar results including tests for different cos- 
mologies. By starting when D(zi n )/D(z = 1) = 0.01 our results 
are converged to better than 1 % for all < z < 1 . 

6. RESOLUTION TESTS 

In order to ensure that our results are properly converged for 
k < 1 ZiMpc -1 between z = 1 and z = we need to understand 
the impact of box size, particle loading, force softening, and 
particle sampling on the numerically determined power spectra. 

6.1. Box Size 

The choice of the box size depends on several factors. In 
principle, one should choose as large a volume as practicable, to 
ensure that the largest-scale modes are (accurately) linear at the 
redshift of interest (in our case between z = 1 and z = 0), improve 
the statistical sampling (especially for BAO), and to obtain ac- 
curate tidal forces. Practical considerations, however, add two 
restrictions to the box size arising from (i) the necessarily finite 
number of particles, and for the PM simulations, (ii) limitations 
on the force resolution. The storage requirements and run time 
for the N-body codes scale (close to) linearly with particle num- 
ber, so running many smaller boxes "costs" as much as running 
one very large box with more particles. However the ability to 
move jobs through the queue efficiently and post-process the 
data all argue in favor of more smaller jobs than one very large 
job. 

The CDM power spectrum peaks roughly at k ~ 0.01 hMpc , 
determined by the horizon scale at the epoch of matter-radiation 
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equality. As the power falls relatively steeply below this value 
of k, a box size of 1 (/z _1 Gpc) 3 , corresponding to a fundamental 
mode of k ~ 0.006 /zMpc -1 , is a reasonable candidate for com- 
paring with linear theory on the largest scales probed in the box. 
[These considerations are of course redshift and erg -dependent: 
at z = 0, small nonlinear mode-coupling effects can be seen be- 
low k ~ 0.1/zMpc -1 (cf. Figure IT). At higher redshifts, these 
effects move to higher k.] Of course, bigger boxes are even bet- 
ter (especially for improved statistics, although this is unrelated 
to linear theory considerations), and a convergence test in box 
size is described below. 

The particle loading is particularly significant as it sets the 
maximum wavenumber below which the power spectrum can 
be accurately determined. As discussed in Section 16*721 the ac- 
curacy of the power spectrum degrades strongly beyond the 
Nyquist wavenumber, which depends on both the box size 
and particle number [see Eqn. ®]. Therefore, a compro- 
mise has to be found between box size and particle loading. 
After having decided the size of the smallest scale of inter- 
est and the maximum number of particles that can be run, 
the box size is basically fixed. In our case, the optimal so- 
lution (considering computational resources) appears to be a 
box size of roughly 1 /r'Gpc on a side and a particle load- 
ing of one billion particles - covering a wavenumber range 
0.0067 /!Mpc"'< k < 3.4 /iMpc" 1 with the upper limit given 
by the Nyquist wavenumber. 

The force resolution for PM codes is a direct function of the 
box size, once the size of the density (or PM) grid is fixed. 
While other codes do not have this restriction in principle, PM 
codes are very fast, and have predictable error properties. In or- 
der to obtain sufficient statistics and accuracy for determining 
P(k), results from many large volume runs at modest resolu- 
tion can be "glued" to those from fewer high resolution runs, 
providing an optimal way to sample the quasilinear and non- 
linear regimes. PM simulations are very well suited to han- 
dling the quasilinear regime; for a Gpc 3 box, a 2048 3 grid pro- 
vides enough resolution to match the high-resolution runs out 
to k ~ 0.5 /zMpc" 1 . 
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FIG. 7. — Comparison of power spectra from two different box sizes, 
1728/r'Mpc (blue: z. = 0, green: z = 1) and 936/i~'Mpc (red: z. = 0, orange: 
z = 1). We show the ratio of the power spectrum at z = and z = 1 to the initial 
conditions at z;„ = 211. All power spectra have been scaled via the growth 
factor to z = 0. For the larger box we average over four realizations, for the 
smaller box over eight realizations so there is more scatter in the small box 
data. The overall agreement of A 2 (k) from the two box sizes is better than 1% 
on scales below k ~ 0. 1 //Mpc~' (the 1% limit is shown by the dotted lines). 
The difference at higher k is due to the difference in resolution of the codes 
which have fixed dynamic range but different box sizes. 



In order to ensure that a Gpc 3 box is sufficient to obtain ac- 
curate results on very large scales, we compare the results from 
several runs in a (936 /r'Mpc) 3 box to a (1728 /z _1 Mpc) 3 box at 
Z = 1 and z = (see Figure|7]i. We divide each realization by its 
initial power spectrum scaled by the growth function to be able 
to compare the two boxes on the largest scales. Each simula- 
tion was run with 1024 3 particles on a 2048 3 grid with one of 
our PM codes. We average results from eight small boxes and 
four big boxes. The agreement between the two boxes is much 
better than one percent. The agreement with linear theory on 
scales below k ~ 0.1 /zMpc -1 is roughly at the percent level and 
much better than this for k ~ 0.01 /zMpc" 1 . We note that for the 
cosmology used in our study, we do not observe a suppression 
of the power spectrum with respect to linear theory by ~ 5% on 
scale s of 0.05 /zMpc" 1 < k < 0.075 /iMpc" 1 as was reported in, 
e.g JSmithetall (120071). 



In a recent paper, Takahashi et al. (2008) discuss finite vol- 
ume effects in detail and propose a way to use perturbation the- 
ory to eliminate these effects. They have two concerns: (i) A 
small simulation volume will lead to enhanced statistical scatter 
on large scales, if only a few realizations are considered, (ii) If 
the simulation volume is too small and the linear regime is not 
captured accurately, the result for the power spectrum will be 
biased low. We overcome the first difficulty by running many 
realizations of our cosmological model. In combination with 
our large simulation volume, we are able to keep the statistical 
noise below the percent level. The second concern is clearly 
valid if the simulation box is too small. With the Gpc 3 and 
larger volumes we consider, no size-related bias is observed. 
The two different box sizes we investigate are in good agree- 
ment as can be seen in F igure [7] One concern with respect to 
the lTakahashi et all (120081) results is that they start their simula- 
tions rather late (zj n = 30) and investigate the results starting at 
Z = 3. As demonstrated in Figure [6] such a late start suppresses 
the power spectrum at quasilinear and nonlinear scales. 

6.2. Mass Resolution 

We investigate the influence of the particle loading on the 
accuracy of the power spectrum by first asking the following 
question: How many particles are required to sufficiently sam- 
ple the density field when calculating the power spectrum? To 
answer this question we start from one of the GADGET-2 simu- 
lations run with a (936 /r'Mpc) 3 box and with 1024 3 particles. 
We determine the power spectrum from this run at z = 0. Next, 
we downsample the 1024 3 particles to 512 3 , 256 3 , and 128 3 
particles by taking the particles which belong to every second 
(fourth, eight) grid point in each dimension. Since the parti- 
cles are downsampled from a fully evolved simulation, evolu- 
tion and sampling issues are separated. 

In the upper panel of Figure[8]the resulting power spectra are 
shown. The lower panel shows the ratio of the power spectra 
from the downsampled distributions with respect to the 1024 3 
particle distribution. In addition, we have marked the Nyquist 
wavenumber divided by two for each power spectrum. The 
Nyquist wavenumber is set by the inter-particle separation on 
the initial grid: 



^Ny - -T- : 



ttN„ 



(9) 



with A p being the inter-particle spacing, N p the cube-root of the 
number of particles, and L, the box size (936 /r'Mpc). Values 
of Jt Ny for the 1024 3 , 512 3 , 256 3 , and 128 3 particle cases are 3.4, 
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FIG. 8. — Importance of particle sampling for calculating the power spec- 
trum. The underlying simulation is the GADGET-2 run with 1024 3 particles. 
All power spectra are measured on a 2048 3 grid. Upper panel: black - power 
spectrum from 1024 3 particles at z = 0, red - from the 512 3 downsampled 
distribution, green - from the 256 3 downsampled distribution, blue - from 
the 128 3 downsampled distribution. Vertical lines denote k-^ y /2 for the three 
cases: 128 3 (black), 256 3 (blue), 512 3 (red). Lower plot: ratios of the down- 
sampled power spectra with respect to the 1024 3 particle power spectrum. The 
dotted line represents the 1 % deviation limit. 
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FIG. 9. — Mass resolution test. The 1024 3 particle ICs have been downsam- 
pled to 512 3 and 256 3 particles and run to z = 0. The plot shows the ratio of 
the obtained power spectra at z = 1 (upper panel) and z = (lower panel) with 
respect to the full 1024 3 power spectrum at z. = (blue, red) and the power 
spectra which are obtained by downsampling the particles at z = (turquoise, 
orange; see also Figure[8j. The vertical lines mark k = &N y /2 for the two cases. 



1.71, 0.86, and 0.43 /zMpc -1 , respectively. As shown in Fig- 
ure[8] all power spectra agree to better than 1% for k < k^ y /2. 
The undersampled particle distributions lead to an overpredic- 
tion of the power spectrum beyond this point due to the increase 
in particle shot noise. As mentioned earlier, a simple shot noise 
subtraction assuming Poisson noise as given in Eqn. © does 
not compensate for this increase. Detailed tests show that the 
shot noise which leads to the overprediction is scale-dependent 
and smaller than Poisson shot noise on the scales of interest. 



(A naive Poisson shot noise subtraction would alter the power 
spectrum at k = 1 /iMpc" 1 at 0.2% at z = and at 1% at z = 1.) 
Thus we are led to conclude that, in the absence of shot noise 
modeling (a difficult and potentially uncontrolled procedure), 
the one-percent accuracy requirement on the power spectrum 
can only be satisfied for wavenumbers, k < k^ y /2. 

The next step is to investigate how the error from an "un- 
dersampled" initial particle distribution propagates through the 
numerical evolution. For this test we first downsample the ini- 
tial particle distribution in the same way as before, at z; n = 21 1, 
from the original 1024 3 particles to 512 3 particles and 256 3 par- 
ticles. We then run the simulations to z = with the same set- 
tings in GADGET-2 as were used for the full run (2048 3 PM grid 
and a softening length of 50kpc). We do not use the 128 3 par- 
ticle set for this test since the corresponding sampling error is 
too large. Results are shown in Figure [9] for outputs at z = 1 
and z = 0. Ratios of the power spectra from the downsampled 
initial conditions (ICs) are shown with respect to: (i) the power 
spectrum from the full 1024 3 run, and (ii) the power spectra 
correspondingly downsampled at z = 1 and z = as shown in 
Figure [8] 

There are two points to note here. First, restricting attention 
to case (i) above, there is a noticeable loss of power below k^ y , 
and second, a steep rise beyond this point. The loss of power 
is not due to the downsampling in the initial condition - as can 
be easily checked by comparing the power spectrum from the 
particles after the IC generation against the desired input power 
spectrum for the given realization - but is due to a discreteness 
effect: a reduction in the linear growth factor from its contin- 
uum value as k — > k^ y . As the evolution proceeds, this suppres- 
sion is reduced due to the addition of nonlinear power, as can 
be seen by comparing the z = 1 and z = results in Figure [9] 
and also by noting the smaller suppression for the case with 
512 3 particles for which the larger k^ y means an enhancement 
in nonlinearity (cf. Figure [8]). The steep rise is a manifestation 
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10 



Precision Determination of the Nonlinear Matter Power Spectrum 



of particle shot noise as can be seen by looking at the results for 
case (ii). For wavenumbers up to k^ y /2 there is no difference 
between the two ratios [case (i) vs. case (ii)] but beyond that 
point the results from case (ii) show a marked reduction (z = 1) 
to almost a removal (z = 0) of the enhancement, consistent with 
the stated hypothesis. We would like to re-emphasize that our 
convergence tests show that a Poisson shot noise subtraction 
alters the power spectrum in the wrong way at the scales of in- 
terest. It enhances the suppression of the power spectrum near 
the Nyquist wavenumber and overcorrects the power spectrum 
at higher wavenumbers. 

The problem we now face is that the (IC downsampling) error 
at k ~ kuy/2 is large: for the 256 3 particle run at z = 1 it is 
~ 20%, and for 512 3 particles it is still ~ 7%. At z = 0, the 
error is ~ 10% for the 256 3 run and ^3% for the 512 3 run. 
Thus, one may wonder if the fiducial 1024 3 particle run can 
itself yield results at k = 1 /iMpc -1 accurate to 1%. 

A brute force approach would be to run with 2048 3 particles 
and check convergence with respect to that simulation. To avoid 
the computational cost of the brute force approach, we take a 
different tack: We extrapolate from the two low-mass resolu- 
tion runs to try and predict the results of the high-mass resolu- 
tion run. The success of Richardson extrapolation when applied 
to power spectra from different force re solution runs has been 
demonstrated by iHeitmann et al] (120051) . We now carry out a 
similar procedure, allowing for both linear or quadratic conver- 
gence. (Details can be found in Appendix iBl) 

Figure [10] shows the results for the extrapolation tests for 
z = 1 and z = 0. Following Eqns. dB5l ) and ( |B7| i, we assume 
linear and quadratic convergence respectively, and predict the 
power spectrum for the 1024 3 particle run, displaying the ra- 
tio of the prediction with respect to the full 1024 3 run. The 
quadratic extrapolation scheme works much better than the lin- 
ear one - out to k ~ 0.8 /r'Mpc the prediction is accurate to 
better than 1%. Obviously, the prediction will not work very 
well beyond the scale set by the mass resolution of the 256 3 



simulation. Nevertheless, the test shows that at k = 1/iMpc 1 
(which is close to k^ y /2 from the 512 3 particle run and below 
£ Ny /2 for the 1024 3 particle run), we will obtain a reasonably 
accurate prediction for a 2048 3 particle run. 

The result for the 2048 3 prediction is shown in FigureQl] At 
z = the 1024 3 particle run has converged at k = 1 /iMpc" 1 to 
better than 1 %, at z = 1 the convergence is slightly worse and 
closer to 2-3%, but here the extrapolation scheme itself is being 
stretched to its limit - the actual result is likely to be better. 
We conclude that our mass resolution will allow a 1% accurate 
calculation at the scale of interest. 

6.3. Force Resolution 

As discussed in Section [3] we employ two A^-body methods 
in this paper: PM simulations with grid sizes of 1024 3 and 
2048 3 and tree-PM simulations. The force resolution of the 
PM runs is insufficient to resolve the power spectrum out to 
k ~ 1 /iMpc -1 (see, e.g., Figure[T4lfor the shortfall of power in 
the PM runs). We therefore discuss only the convergence prop- 
erties of the tree-PM algorithm out to k ~ 1 h Mpc -1 . Since the 
GADGET-2 runs with 1024 3 particles are computationally ex- 
pensive, and the force softening primarily affects small scales, 
we chose to downscale the simulation box and number of par- 
ticles for this test to 256 3 particles in a 234 /T 1 Mpc box (a re- 
duction by a factor of 64 from the main runs). Following the 
practice in the larger runs, the PM force grid is set to twice 
the number of particles in one dimension, resulting in a 512 3 
PM mesh. All the other code settings are the same as for the 
large runs and we vary only the force softening to test for the 
effects of finite force resolution. The effective force resolution 
lengths range from 400 kpc to 25kpc (50kpc is used in the large 
runs). The results for z = and z = 1 are shown in Figure[T2] At 
k ~ 1 /zMpc -1 , the difference between 50kpc and 25 kpc is well 
below 0.1% for both redshifts, and therefore comfortably within 
our requirements. In fact, meeting the force resolution require- 
ments at k ~ 1 h Mpc -1 with the tree-PM algorithm is computa- 
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FIG. 11. — Mass resolution convergence study with respect to the extrapo- 
lation results at z = 1 and z = 0. The content of this figure is very similar to 
that of Figure [To] but here we show results with respect to the 2048 3 particle 
prediction from the (quadratic) Richardson extrapolation. This allows us to 
investigate the convergence properties of the 1024 3 particle runs (see text for 
details). The vertical lines are the same as in Figure|9] 
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FIG. 12. — Force resolution convergence study at z = 1 and z = with 
GADGET-2. The 512 3 PM grid is the same in all five runs, and the force 
resolution is varied between 25 kpc and 400 kpc. At k ~ 1 AMpc -1 , a force 
resolution of 100 kpc already leads to results converged well below 1% at both 
redshifts with respect to the 25 kpc resolution run. 
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tionally much less demanding than meeting the mass resolution 
requirements. It may be that for power spectrum simulations a 
hybrid or adaptive PM code is the most computationally effi- 
cient route, though other uses of the simulations may be more 
sensitive to resolution. 

The size of the PM mesh is a separate issue, and significant 
in its own right. If high accuracy is desired the mesh should not 
be chosen to be too small, as this increases the PM error and 
pushes the handover between the tree and the mesh to larger 
scales. In tests carried out to determine the size of the PM grid, 
we observed an unphysical suppression of the early-time power 
spectrum at quasi-linear scales for the smaller meshes. 

6.4. Time Stepping 

Most iV-body codes use low-order - typically, second or- 
der - symplectic time-stepping schemes. (Full symplecticity 
is not achieved when adaptive time-stepping is employed.) The 
choice of the time variable itself can vary, although typically it 
is some function of the scale factor a, e.g., a itself or the natural 
logarithm of a. PM codes most often use constant time stepping 
in a or In a. Higher-resolution codes use adaptive, as well as in- 
dividual particle time-stepping. Hybrid codes that mix grid and 
particle forces such as tree-PM, have different criteria for time- 
stepping the long-range forces as compared to the short-range 
forces, where individual particle time-steps are often used. Be- 
cause of these complexities, it is important to check that the 
time-stepping errors are sub-dominant at the length-scales of 
interest for computing the mass power spectrum. 

The GADGET-2 runs in this paper use In a as the time vari- 
able. The PM calculations within GADGET-2 use a global time 
step; we found 256 time steps sufficient for this part. The tree 
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FIG. 13. — Time stepper convergence for A 2 (k) using lineal' (upper plot) and 
logarithmic (lower plot) time-stepping in a at z = 0, as a function of number 
of time steps. The k value chosen is for the largest mode in the box, k = 
6.7 X 10~ 3 AMpc -1 (black stars). The black triangle shows the result from a 
GADGET-2 run with adaptive time stepping in In a, the blue box is the power 
spectrum from the initial condition scaled by the linear growth factor to z = 0, 
the red circle the A 2 (k) value for time-stepping linear in a extrapolated to zero 
assuming quadratic convergence, and the turquoise cross the same quantity for 
In a. All (extrapolated) values from the simulations agree with linear theory 
to 0.2% or better, the simulations themselves agreeing to better than 0.04% 
taking the GADGET-2 run as the reference. The pink line shows a quadratic 
fit to the data points. 



algorithm for the short-range forces uses an adaptive time step- 
ping scheme and our runs use a total of about 3000 time steps. 
The criterion for the adaptive time stepping is coupled to the 
softening length e via : Af = y/2rje/\a\ where r\ allows adjust- 
ments in the time stepping; we use r\=\% (note that here a is 
the acceleration). 

The basic convergence test is to consider the largest mode in 
the box, in this case, k = 6.7 x 10~ 3 /iMpc _l and compare the 
numerical results for P(k) with that expected from linear the- 
ory, which should be very accurate at these very large scales. 
We chose to investigate both time variable choices, In a and a. 
The results are shown in Figure Q~3] All the test runs are in 
pure PM mode on a 1024 3 grid, with the tree switched off in 
GADGET-2 (there is no need for high force resolution in this 
test). For time-stepping linear in a we show results for roughly 
600 and 1,200 time steps, for the time stepper in In a we show 
results for Alna » 0.005, 0.01, 0.02, 0.04, and 0.08. In addi- 
tion, we fit two curves through the results assuming linear and 
quadratic convergence. As expected from a second order in- 
tegrator, the quadratic fit is in very good agreement with the 
data points. Quadratic extrapolation of the results for the two 
time stepping schemes from finite k to zero is in very good 
agreement with linear theory, to better than 0.2% - about the 
deviation expected given the dimensionless power at the fun- 
damental mode of the box. If we take the adaptive time step 
run as the reference (rather than linear theory), the agreement is 
better than 0.04%. Adaptive time-stepping is expected to yield 
results very close to In a stepping on large scales, since for the 
long-range force even the adaptive time-stepper run is constant 
in lna with Alna = 0.02. The excellent agreement with time- 
stepping in a confirms the robustness of the different schemes. 
Since our interest is in generating the power spectrum at percent 
accuracy at minimal computing cost, we conclude that the lna 
time-stepping scheme with approximately 250 time steps is a 
good compromise for the PM runs to obtain an accurate power 
spectrum at large scales. 

7. MATCHING LOW AND HIGH RESOLUTION POWER SPECTRA 
AND COMPARISON WITH HALOFIT 

Last, we comp are our results wi th the standard fitting for- 
mula, HaloFit dSmith et aI1l2003l). currently used for analy- 
sis of e.g. weak lensing dat a dJarvis et al. 20061: iMassev et al.l 
120071; iBeniamin et ail 120071; iFu et all 120081) or for forecasts 
on the impro vement of c osmological constraints from future 
surveys dTang et alj|2008l) . HaloFit provides the nonlinear 
power spectrum over a range of cosmologies in a semi-analytic 
form. It is based on a combination of t he halo model approach 
(for a review of the halo model, see e.g.. lCoorav &Sheth2002) 
and an an alytic description of th e evolution of clustering pro- 
posed by Hamilton et al.l d 199 lb . In addition, the fit is tuned 
to simulations by introducing two new parameters: an effective 
spectral index on nonlinear scales, n e ff, and a spectral curvature 
C. The combination of analytic arguments and tuning to results 
from iV-body simulations has led to the most accurate fit for the 
nonlinear power spectrum to date (as we will show below, the 
fit is accurate to ~ 5- 10%). As mentioned above, we use here 
the CAMB implementation of HaloFit. 

In order to compare simulation results to a smooth fit, we 
first combine 16 realizations from the PM runs in the low k re- 
gion with one high-resolution run, as shown in Figure [14] At 
around k = 0.6 /iMpc" 1 the lower resolution of the PM runs be- 
gins to become apparent and the result falls below that from 
GADGET-2. Conservatively, we match the two power spectra at 
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FIG. 14. — Matching of an ensemble of low resolution runs with one re- 
alization of a high-resolution GADGET-2 run. The upper panel shows the 
average from 16 realizations from the low resolution PM runs (red) and the 
power spectrum from the GADGET-2 run (black). The lower panel shows 
the ratio of the low resolution ensemble with respect to the GADGET-2 run. 
Out to k ~ 0.5 /iMpc~' the difference is less than one percent (disregarding 
the noise from the single realization). We match the two power spectra at 
k ~ 0.3 /iMpc~' , at which point the noise in the single realization is small 
enough, yet the resolution of the PM runs is sufficient to accurately resolve the 
power spectrum. 
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FIG. 15. — Ratio of two realizations at initial and final redshift. Both sim- 
ulations are carried out with GADGET-2 at the standard setting. Results as 
shown have been smoothed by averaging over every five ^-values. The ratio 
beyond our matching point of low and high resolution simulation is at the per- 
cent level, confirming that one realization of a high-resolution run is sufficient. 



k = 0.3 /zMpc" 1 . At this point, the variance from the single real- 
ization of the GADGET-2 run is small enough that the matching 
leads to a smooth power spectrum. One concern might be that a 
single realization is insufficient to capture the behavior on small 
scales accurately: Because of mode coupling it is not obvious 
that fluctuations on large scales do not also cause substantial 
effects on small scales. In Figure [15] we show that, due to the 
large box size, this is not a concern at least at the percent level 
of accuracy. The figure shows the ratio of two different realiza- 
tions at the initial and final redshift. Both simulations are run 
with GADGET-2 at our standard settings. The variations beyond 
k ~ 0.3 /zMpc" 1 are at the percent level and appear to be free of 
systematic trends. 

The ratio of the matched power spectrum to the predic- 
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FIG. 16. — Comparison of the simulation power spectrum to HALOFlT. 
Shown is the ratio of HALOFlT with respect to the simulation result. The simu- 
lation result has been obtained by combining the PM runs and the GADGET- 
2 ran at k = 0.3 /iMpc -1 and it has been smoothed by averaging over every 
five ^-values to reduce the noise for the comparison. The HALOFlT result is 
approximately 5% lower than the result from simulations. 



tion from HaloFit is shown in Figure [16] In this case, the 
HALOFlT prediction falls roughly 5% below the simulation. 
The procedure for combining the simulation results can be seen 
to work very well, as there is no discontinuity at k = 0.3 /iMpc -1 
from the matching . Our result is in good agreement with, e.g., 
ISmith et al.l d2008h as well as Ma (2007), who find a 5% supres- 
sion for HaloFit at k ~ 0.1 /zMpc -1 . At larger k, however, the 
results in Ma (2007) may not be very accurate, due to limita- 
tions in force resolution. 

8. CONCLUSION AND OUTLOOK 

The advent of precision cosmological observations poses a 
major challenge to computational cosmology. With observa- 
tional results accurate to the percent level a significant uncer- 
tainty in extracting cosmological information from the data 
is due to inaccuracies in theoretical templates. At the re- 
quired level of accuracy large scale simulations are unavoid- 
able, since the nonlinear nature of the problem makes it impos- 
sible to derive analytic or semi-analytic expressions for statis- 
tics such as the matter power spectrum, at an accuracy better 
than ~ 10%. While simulations in principle should yield results 
at sub-percent accuracy, in practice this is a non-trivial task due 
to uncertainties in the numerical implementation and modeling 
of relevant physical processes. 

In this paper, we showed that it is possible to obtain the non- 
linear matter power spectrum at sub-percent/percent accuracy 
out to k ~ 1 /zMpc -1 between z = 1 and z = 0. This wavelength 
regime is important for ongoing and near-future weak-lensing 
surveys. The restriction to these (large) length scales has two 
major advantages: baryo nic eff e cts are subdominant on thes e 
scales as shown by, e g.. IWhitd (120041) : IZhan & Knox] d2004 : 
Uing et al.1 (120061) and iRudd et"atl (120081) . and the numerical re- 
quirements in this regime remain rather modest. Each simula- 
tion can be carried out in a matter of days on parallel comput- 
ers with several hundred processors. Pushing the scales beyond 
k^\h Mpc -1 will require advances in our understanding of the 
implementation of baryonic physics as well as advances in al- 
gorithms and computational power. 

Baryon acoustic oscillations, on the other hand, need predic- 
tions on larger scales, but at even higher accuracy levels than 
we have targeted in this paper. For BAO, predictions approach- 
ing the 0.1% level are required. Our current work should be 
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helpful in guiding such an effort, though the simulation volume 
will have to be increased roughly by an order of magnitude to 
suppress finite sampling variance. 

The aim in this paper was to derive a set of numerical re- 
quirements to obtain an accurate power spectrum by performing 
a large suite of convergence and comparison tests. As shown 
here, the simulation volume and, especially, the particle loading 
are two major concerns in obtaining an accurate power spec- 
trum. The simulation volume has to be in the ^Gpc 3 range, 
leading to a minimum requirement of ~ 1 billion particles. Fur- 
ther increase in volume would be helpful, but would require a 
concomitant increase in the number of particles, greatly adding 
to the computational burden. The 1 Gpc 3 /1 billion particle sim- 
ulation is a good compromise between sufficient accuracy and 
computational cost. 

Besides a large simulation volume and good particle sam- 
pling, initialization of the simulation also plays an important 
role. To guarantee converged results, the simulation must be 
started at a high enough redshift. We found that a starting red- 
shift of z; n = 200 is sufficient to get accurate results between 
Z - 1 and z = 0. 

We found that results for the power spectrum are rather stable 
to changes in the number of time steps. This is clearly related to 
the fact that our resolution demands are relatively modest. For 
the PM runs, a few hundred time steps are sufficient, while for 
the tree-PM runs the overall number of time steps is a factor of 
ten larger. We emphasize that the simulation settings discussed 
here will lead to the required accuracy only up to k ~ 1 /iMpc -1 . 
While these settings can be used as a guideline for other simu- 
lation aims, they do not replace convergence tests that must be 
performed for each new problem, if one desires high precision 
results. 

Having established the ability to generate power spectra with 
sufficient accuracy from A^-body simulations, the next major 
question that arises is how to use these costly simulations for 
parameter estimation, e.g., via Markov Chain Monte Carlo. 
To address this problem, w e have recently introduced the cos- 
mic calibration framework dHeitmann et alj[2006h lHabib et"al] 
l2007t ISchneideretaT]|2 008 ) which is based on an interpolation 
scheme for the power spectrum (or any other statistic of inter- 
est) derived from a relatively small number of training runs. 

The next step in generating precise predictions for the mat- 
ter power spectrum is to determine the minimum number of 
cosmological models needed to build an accurate emulator and 
then to construct the emulator from a set of high-precision sim- 
ulations. In the second paper of this series we establish that 
30-40 cosmological models are sufficient to explore the param- 
eter space for wCDM cosmologies (constant w) given the cur- 
rent constraints on parameter values. The third and final paper 
will present results from the simulation suite designed and dis- 
cussed in the second paper, and will include a power spectrum 
emulator that will be publicly released. 
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APPENDIX 

ZEL'DOVICH APPROXIMATION VERSUS 2LPT 

The initial conditions for A^-body simulations are usually 
generated by displacing particles from a regular grid using the 
Zel'dovich approximation. This amounts to a first order expan- 
sion in Lagrangian perturbation theory. The use of a higher 
order Lagrangian approximation scheme to set up init i al con - 
ditions has been suggested recently, e.g., ICrocce et alj (|2006). 
Such a scheme is formally higher order in accuracy, compared 
to the Zel'dovich approximation, and has transients which de- 
cay much faster with the expansion of the Universe (a~ 1 rather 
than a" 1 ). In the 2LPT formalism, the particle displacement is 
obtained in second order Lagrangian perturbation theory, an ad- 
ditional contribution is added to the Zel'dovich approximation 
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FIG. A17. — Upper panel: distribution of the initial displacements of all par- 
ticles at different starting redshifts (z: m = 200, 100, 50, 25). The displacement 
is measured with respect to the mean inter-particle spacing. For z m = 200, the 
rms displacement is approximately 0.05, while for Zi a = 25 it increases by a 
factor of ten. Lower plot: 2LPT correction. The distributions show the ad- 
ditional contribution in the initial move to the Zel'dovich approximation. For 
Zm = 200 this additional move is on average 4 ■ 10~ 5 and for z in = 25 it is 0.004 
of the mean inter-particle spacing. In both cases this is a small fraction with 
respect to the Zel'dovich move. In both plots, the y-axis is scaled with respect 
to all particles. 
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FIG. Al 8. — Power spectra ratios for four different initial redshifts. The ini- 
tial power spectrum obtained from Zel'dovich initial conditions is divided by 
the power spectrum from the 2LPT initial conditions. Overall, the Zel'dovich 
initial conditions have slightly less power on the smallest scales. Results are 
shown out to &Ny/2 = irNp/lL = 1,57/jMpc . Remarkably, even if the initial 
conditions are generated as late as z m = 25, the difference in the power spectra 
is below 1 % at the smallest scales. For zin = 200, the difference on all scales is 
far below one percent. 



given in Eqn. © leading to: 

x(q) = q-D 1 V ?( /. (1) +Z) 2 V, 
c/x 



AC*) 



= = -D x f x HV q 4> {l) +D 2 f 2 HV q ^ 2 \ 



where </> (2) is obtained from solving 
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and D 2 is the second order growth function. In the following, 
we investigate the contributions from the second terms in the 
p ositions and v e locit ies of the particles at different redshifts. 

ICrocce et al.l d2006f) have made a serial 2LPT code publicly 
available. Their code uses approximations for the growth func- 
tions in first and second order. (By contrast, the initialization 
routine used for this paper solves the differential equation for 
the linear growth function directly, without making approxima- 
tions.) For a ACDM cosmology these approximations are given 
by: 
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with t bein g conformal t i me. T he approximation for D x can 
be found in ICarroll etal] d!992h . For f x and f 2 the following 
approximations are made: 

/i«^ 9 , / 2 « (A6) 
A detailed discussion of the exact differential equations for the 
growth function up to th ird order and the rel iability of these ap- 
proximations is given in Bou chet et alj d 1995b . In order to limit 
computational expense, we restrict our tests using this code to 
256 3 particles in a 256/r'Mpc volume. This choice is suffi- 
cient to study the general question, as the inter-particle spacing 
is the same as in the main runs. In keeping with our general 
philosophy of redundancy and cross-checking we also indepen- 
dently implemented a 2LPT initial conditi ons generator whic h 
gave essentially the same results as that of lCrocce et alj (2006). 
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FIG. A19. — Ratio of histograms of the three velocity components from 
the Zel'dovich approximation and the 2LPT approach. The insets show the 
regimes between -1000 km/s and 1000 km/s where the large majority of the 
particles reside. Here the difference is sub percent. The different colors rep- 
resent different starting redshifts, the difference becoming smaller for higher 
redshift starts. 



We generate four sets of initial conditions at z; n = 200, 100, 
50 and 25. All of the initial conditions have the same phases 
and can therefore be compared directly. First, we measure the 
displacement from the Zel'dovich approximation; results are 
shown in the upper panel of Figure [ATTI For this one realiza- 
tion, the rms displacement at z Xn = 200, which is the starting 
redshift for our main simulations, is around 5% of the mean in- 
terparicle spacing. By delaying the start until z; n = 25, the rms 
displacement grows by a factor of ten. The 2LPT correction, 
given by the second term in Eqn. (IAU , is negligible at Zm = 200, 
being smaller than 10~ 4 on average. In fact at this point numeri- 
cal accuracy might be questioned, since the approximations for 
the growth functions might not be accurate at this level. Fig- 
ure IA18I shows the ratio of the initial power spectra from the 
Zel'dovich and the 2LPT approximations. As for the displace- 
ments, convergence with increased redshift is very apparent. At 
a starting redshift of z; n = 200, both power spectra agree to bet- 
ter than 0.02%. Even starting at very late times (z; n = 25) only 
leads to a 1 % difference between the initial power spectra. 

Next we measure the differences in the initial velocities from 
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FIG . A20 . — Comparison of the power spectra from simulations started from 
Zel'dovich initial conditions and 2LPT initial conditions at z = 1 (upper panel) 
and z = (lower panel). Shown are the ratios of power spectra from starts at 
redshift 1 +z in = 200, 100, 50, 25. The start at l+z in = 200 leads to an agreement 
of the power spectra better than 0.5% at k ~ 1 AMpc -1 and better than 0.2% 
at z = 0. At larger scales, k < 0. 1 h Mpc~' , the agreement is basically perfect. 
Therefore, the Zel'dovich initialization scheme started at 1 +z m = 200 fulfills 
our accuracy requirements comfortably. 



the two approximations. The results are shown in Figure |AT9l 
We display the three velocity components v x , v y , and v, sepa- 
rately. The main difference occurs in the tails of the velocity 
distributions. Independent of redshift a negligible number of 
particles (fewer than 0.5%) live in these tails with absolute ini- 
tial velocities larger than lOOOkm/s. Ignoring these tails (see 
the insets in Figure lAl9b . the difference in the velocities be- 
tween 2LPT and ZA starting at different redshifts is below 1%. 
At Zi n = 200, the difference is less than 0.1%. At this precision, 
the inaccuracy from the approximations for the growth function 
at first and second order is probably larger than the error from 
the Zel'dovich approximation. 

The velocity differences are highly correlated with density 
however (see also Figure [5]), and to understand this effect we 
evolve initial conditions created from the ZA and 2LPT forward 
to z = 0. We use our parallel 2LPT code, which does not rely 
on an approximation for the growth function, to generate initial 
conditions with 5 12 3 particles in a 468 /r'Mpc box - downscal- 
ing our main runs by a factor of eight. The ICs are generated at 
four different redshifts, 1 +z; n = 200, 100, 50 and 25 and evolved 
to z = using a tree-PM code. We measure the power spectrum 
of the evolved particles at z = 1 and z = 0. The results are shown 
in Figure IA20I where we see a shortfall in power at high k in 
the ZA starts as compared to the 2LPT starts but convergence 
as Zi n is increased. At k ~ 1 /iMpc -1 the evolved power spectra 
from both sets of initial conditions at 1 + z; n = 200 show excel- 
lent agreement, better than 0.5% at z = 1 and 0.25% at z = 0. 
We therefore conclude that our starting redshift, 1 +z m = 200, is 
high enough to avoid any problems arising from possible inad- 
equacies of the Zel'dovich approximation. 

RICHARDSON EXTRAPOLATION 

Richardson extrapolation is a method to compute the limiting 
value of a function that is assumed to have a smooth behavior 
for small deviations around the evaluation point. Suppose we 



have such a function /, then it is plausible to assume that 

/(0+A) = /(0) + c 1 A + c 2 A 2 + c 3 A 3 + ---. (Bl) 

For many quantities derived from numerical simulations, it is 
not often a priori obvious what the convergence structure, i.e., 
the values of the coefficients, c,, happens to be, even to the ex- 
tent of knowing which of the coefficients are zero or non-zero. 
Nevertheless, for small enough values of the deviation, A, one 
can numerically establish the values of the leading order coeffi- 
cients. By doing so, one can then go on to improve estimates for 
the desired limiting value /(0) using Richardson extrapolation. 

As a simple example, consider the case of non-zero c\ (linear 
convergence) for some quantity, say the power spectrum at a 
given value of k, as a function of the mesh spacing in a PM 
code. Then, if we write, for a 256 3 mesh, 



/(2A)~/(0)+2dA, 
for512 3 and 1024 3 meshes we would have, 
/(A)~/(0) + Cl A, 
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where A has been taken to be the mesh spacing for the 512 3 
grid. Eqns. ( IB2b and (IB3b then predict an estimated value for 
the 1024 3 run 



5/(A)-i/(2A) 



(B5) 



which can be used to test whether linear convergence is holding 
for the particular range of values of A. If the test is successful, 
one can then proceed to obtain an estimate for the continuum 
prediction (A = 0) from the 512 3 and the 1024 3 simulations, via 

' A N 
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For the case of quadratic convergence (ci = 0, c 2 ^0), the 
extrapolation from the 256 3 and the 512 3 mesh to the 1024 3 
mesh reads: 
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and the estimate for the continuum from the 512 3 simulation 
and the 1024 3 simulation is given by: 



4 /A 2 \ 1 a2 
/(0) ~ -/ ( — J - -/( A )■ 



(B8) 



Given 3 simulations one can choose to estimate two non-zero 
coefficients, and test the assumed convergence model. 
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